COMMENT ⊗   VALID 00006 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002
C00004 00003
C00008 00004	IF EQU(ASK_STRING(" SKIP FALSE POSITIVES "),"Y") THEN GO TO SECOND_LOOP
C00010 00005	for a range of conf values, find the corresponding 1% and 5%
C00013 00006	SECOND_LOOP:
C00014 ENDMK
C⊗;
comment
	this version has only the test for false positives
	ie finds the 1% and 5% limits for zero signal;

BEGIN "EDGE_TEST"
	comment
		find best choice of conf and diff for Hueckel edge
		operator based on chosen false positive rate for
		no signal.  For a given false positive rate, find
		choice which requires the least height step to
		achieve a 50% false negative rate.  This uses
		steps generated by a Monte Carlo procedure;

DEFINE SAFE_ARRAY=1;	comment this makes all arrays safe;
REQUIRE "VISION.2[SYS,TOB]" SOURCE_FILE; comment vision language package;
REQUIRE EDGE_FILE;
REQUIRE EDGE_PRINT;
REQUIRE TV_VARIABLES;
REQUIRE "STEP4.SAI[FAK,TOB]" SOURCE_FILE;

comment
	conf=sqrt(s0sq/(s0sq+n0sq)) ie cosine
	diff=s0sq;

BOOLEAN SIMPLE PROCEDURE YES(STRING S);
BEGIN "YES"
S←ASK_STRING(S);
IF EQU(S,"Y")∨EQU(S,"y")∨EQU(S,"YES")∨EQU(S,"yes")
THEN RETURN (1) ELSE RETURN (0);
END "YES";
DEFINE NBHD=8;	comment the size of the array with fake intensity values;
LABEL FIRST_LOOP; ! entry for false positives and others;
LABEL SECOND_LOOP; ! entry for false negatives;
LABEL ERR_TEST; ! entry for angular and position errors;
LABEL FALSE_NEG; ! entry for false negatives;
LABEL AVG_DIFF;  ! entry for avg value of diff;
LABEL HISTO_LOOP; ! loop entry for histograms repeat with change scale;
DEFINE MAX_COUNT=2000; ! the max number of test samples generated;
INTEGER_ARRAY INTENSITY[1:NBHD*NBHD]; ! store fake intensities here;
REAL_ARRAY N0SQA[1:MAX_COUNT+20];	! store n0sq from edge oper here;
REAL_ARRAY S0SQA[1:MAX_COUNT+20] ;	! store s0sq from edge oper here;
INTEGER_ARRAY COUNT[0:100];	! histogram of scaled s0sq values;
SHORT INTEGER K,K1;
INTEGER	ENOUGH; ! how many samples (< max_count);
SHORT REAL STEP_MEAN,STEP_AMPLITUDE,NOISE_STD;
	comment the mean,amplitude,noise of fake steps;
SHORT REAL NSQ; ! scaling factor which depends on noise squared;
SHORT REAL NOISE_SCALE;	! scaling factor for histograms;
SHORT REAL THETA,RHO; ! angle,perp dist of edge from origin;
SHORT INTEGER BYTES;	! number of bytes per word;
SHORT REAL LAMBDA; ! wavelength (used to get phase near an edge);
SHORT INTEGER N1,N5,KN1,KN5; comment 1% and 5% level counts and indices;

ENOUGH←MAX_COUNT;

comment now init the edge operator;
CONF←.05;DIFF←.05; ! choose low values to reject nothing;
BITS←6;
MODIFY(BITS);
LINLEN←WORDS(NBHD,BITS);
TVWORD←GIOWD(INTENSITY);
LSIDE←1;RSIDE←NBHD;
FLINE←1;LLINE←NBHD;
EJINIT;
DIFF←0;CONF←.05;

comment set up values for fake edges with zero height;
NOISE_STD←5.0; ! amplitude of gaussian noise;
STEP_MEAN←(2↑BITS-1)*.5; ! offset to get positive intensities;
STEP_AMPLITUDE←0; ! step height for fake edges;
NOISE_SCALE←20.0;

IF EQU(ASK_STRING(" SKIP FALSE POSITIVES "),"Y") THEN GO TO SECOND_LOOP;
FIRST_LOOP:
MODIFY(NOISE_STD);
MODIFY(ENOUGH); ! ask user if he wants another value;
IF ENOUGH>MAX_COUNT THEN ENOUGH←MAX_COUNT;
MODIFY(STEP_AMPLITUDE);

OUTSTR(CRLF&"FALSE POSITIVES ");
FOR K1←1 STEP 20 UNTIL ENOUGH DO BEGIN "ACCUM1"
FOR K←K1 STEP 1 UNTIL K1+20 DO BEGIN "ACCUM" ! report every 20 events;
	SHORT INTEGER N;
	LABEL CALLFK;
CALLFK:	FAKE(NBHD,NBHD,THETA,RHO,LAMBDA,
	STEP_MEAN,STEP_AMPLITUDE,NOISE_STD,NOISE,INTENSITY);
		comment generate an array with zero edge,random;
	EJLI(NBHD/2,NBHD/2,0,0); ! try the edge operator on it;
	N0SQA[K]←N0SQ;	! save the noise squared;
	S0SQA[K]←S0SQ;	! and save the signal squared;
	END "ACCUM";
PRINT(K);
END "ACCUM1";
comment for a range of conf values, find the corresponding 1% and 5%
false positives (above threshold with no edge signal) for fake
edges whose noise sq and signal sq have been saved;
N1←ENOUGH%100;
MODIFY(N1);
TCRLF;PRINT(N1);
HISTO_LOOP:
MODIFY(NOISE_SCALE);
NSQ←NOISE_SCALE/(NOISE_STD↑2);
FOR CONF←.95 STEP -.05 UNTIL .40 DO BEGIN "CONF_LOOP"	! loop over conf;
	SHORT REAL CONFSQ;
	SHORT INTEGER N;
	CONFSQ←CONF↑2;
	comment now form histograms over scaled s0sq values;
	FOR K←0 STEP 1 UNTIL 100 DO COUNT[K]←0; ! zero the histogram;
	FOR K←1 STEP 1 UNTIL ENOUGH DO BEGIN "ACCUM"
	SHORT REAL S;
	S←S0SQA[K];
	IF S/(S+N0SQA[K])>CONFSQ THEN BEGIN ! if above conf;
		N←S*NSQ+.5;
		IF N>100 THEN N←100;
		COUNT[N]←COUNT[N]+1;
		END;
	END "ACCUM";
OUTSTR(CRLF&" OUTPUT HISTOGRAM");
PRINT(CONF);
	INCHWL;
IFC 0 THENC
comment print histogram;
FOR K←0 STEP 1 UNTIL 100 DO BEGIN
	TCRLF;PRINT(K);OUTSTR(" "&CVS(COUNT[K]));
	END;
ENDC;

comment now find the 1% and 5% levels of the histogram tails;
N←0;
KN1←100;
FOR K←100 STEP -1 UNTIL 0 DO BEGIN
	N←N+COUNT[K];
	IF N<N1 THEN KN1←K;
	END;
OUTSTR(CRLF&"DIFF CUTOFF LIMITS "&CVG(N1/ENOUGH));
OUTSTR(CRLF&CVG(KN1/NSQ));
OUTSTR(CRLF&"DIFF CUTOFF LIMITS IN UNITS OF SQUARE STD DEV");
OUTSTR(CRLF&CVG(KN1/NOISE_SCALE));
PRINT(KN1);
END "CONF_LOOP";
IF EQU(ASK_STRING(CRLF&"REPEAT HISTOGRAMS "),"Y") THEN GO TO HISTO_LOOP;
IF EQU(ASK_STRING(CRLF&"REPEAT VALUES"),"Y") THEN GO TO FIRST_LOOP;
SECOND_LOOP:
END "EDGE_TEST";